Bottom-up derivation of an effective thermostat for united atoms simulations of water 
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In this article we derive the effective pairwise interactions in a Langevin type united atoms model of water. 
The interactions are determined from the trajectories of a detailed molecular dynamics simulation of simple 
point charge water. A standard method is used for estimating the conservative interaction, whereas a new 
"bottom-up" method is used to determine the effective dissipative and stochastic interactions. We demonstrate 
that, when compared to the standard united atoms model, the transport properties of the coarse-grained model is 
significantly improved by the introduction of the derived dissipative and stochastic interactions. The results are 
compared to a previous study, where a "top-down" approach was used to obtain transport properties consistent 
with those of the simple point charge water model. 
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I. INTRODUCTION 

Water is the most important biological solvent in natural 
systems. Accurate molecular dynamics modeling of, for ex- 
ample proteins or membranes, depend critically on how well 
the model captures the properties of the ambient water. Un- 
fortunately, the structure of the water molecule, especially its 
tendency to build local hydrogen networks, makes it hard to 
accurately model liquid water. A large number of different 
molecular dynamics (MD) models have been proposed, each 
with advantages and shortcomings, and often tailored to fit a 
particular use. See e.g. 01 |H for an overview. In addition 
to these difficulties, all atomistic models of water share the 
problem of being too computationally demanding for model- 
ing the time and length scales that are relevant for many self- 
assembly and folding processes in biochemistry. To circum- 
vent these limitations several coarse-grained models of water 
have been suggested in the literature, for example the united 
atoms model that we are considering in this paper. While it is 
clear that many central properties of water can only be mod- 
eled implicitly in the aggregated models, especially properties 
depending on the extended hydrogen networks, the models' 
computational efficiency enables the study of processes which 
occur on scales unattainable using more detailed models. 

In this article we consider how to determine the effective 
deterministic and stochastic interaction potentials in a united 
atoms representation of water. Generically, a coarse-grained 
model of a mechanical system is described by Langevin dy- 
namics. Under certain conditions, for example when the de- 
grees of freedom removed by the coarse-graining are aver- 
aging, the resulting dynamics is deterministic, but in general 
both noise and dissipation must be included in the reduced 
representation. It is somewhat surprising that this fact is of- 
ten ignored when effective potentials for united atoms models 
of molecular systems are derived. Under the assumption that 
the effective interactions are pairwise, an important result of 
Henderson guarantees a one-to-one correspondence between 
the pair correlation function and the pairwise potential 01 . 



Using this result, measurements of center of mass pair corre- 
lations from detailed molecular dynamics simulations can be 
used to determine the conservative potential to be used in a 
united atoms simulation Jl II II H, 11 El (alternatively, di- 
rect time averaging over the fast degrees of freedom can be 
used, see e.g. Refs. ITllfl2l) . Importantly, the equilibrium pair 
distribution is not altered by the noise and dissipative interac- 
tions foil . This is perhaps the reason why these interactions 
are often ignored or represented by a weak general purpose 
thermostat. 

The elegant result by Henderson also causes problems by 
the simple fact that it is completely general. Since the pair cor- 
relation function can usually be measured in a detailed simu- 
lation, we can also usually derive an effective pairwise po- 
tential, even if the coarse-grained representation is not closed 
and therefore has no meaning. The conclusion is that requir- 
ing the coarse-grained interactions to recreate the correct pair 
distribution does not guarantee any correspondence between 
the dynamic properties of the detailed system and the coarse- 
grained system. We may even argue that it is not clear if 
it is worth while to simulate the exact time dynamics of the 
coarse-grained model; If the correspondence between the de- 
tailed and the coarse-grained system is only their respective 
equilibrium configurations, it is probably more efficient to use 
a Monte Carlo algorithm to minimize the free energy and to 
sample from the ensemble of equilibrium configurations. 

We have demonstrated in earlier work that the transport 
properties of simple point charge water (SPC), simulated with 
MD, and the corresponding united atoms model, based on 
Henderson's theorem, differ significantly Hill . In the same 
study it is was shown that by adding a DPD thermostat to the 
united atoms model, the transport properties can be tuned to 
correctly represent those of the detailed MD simulation, while 
at the same time preserving the equilibrium pair correlation 
function. The motivation for the DPD thermostat is that it is 
the most general pairwise interaction that respects the local 
symmetries in a mechanical system, i.e., conservation of mo- 
mentum and angular momentum. The disadvantage of this ap- 
proach is clearly that the correspondence between the detailed 
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dynamics and the coarse-grained dynamics may be lost. 

In a recent publication we introduced a method for deter- 
mining the stochastic and dissipative interactions in systems 
simulated with a dissipative particle dynamics (DPD) ther- 
mostat [13]. The method can be viewed as a compliment to 
the inverse Monte Carlo technique JH 0] used to reconstruct 
the conservative interaction based on Henderson's result. To- 
gether, these two methods can be used to determine the de- 
tailed structure of the effective conservative, dissipative, and 
stochastic interactions in a united atoms model. In addition, 
the framework gives a clear indication on whether or not the 
united atoms representation is a valid coarse-graining. 

In this paper, we apply the above framework to a united 
atoms model of SPC. This is possible using the assumption 
that the reduced degrees of freedom are not averaged, but 
rather leads to noise and dissipation that can be represented by 
the stochastic and dissipative interactions in a DPD thermo- 
stat. This assumption should be viewed as similar to the pair- 
wise potential assumption in Henderson's theorem. The trans- 
port properties of the resulting united atoms-DPD model are 
much improved compared to the original united atoms model. 
We discuss the implications of this result for the understand- 
ing of united atoms models and how well a DPD thermostat 
manages to represent the reduced degrees of freedom. 



II. THEORETICAL BACKGROUND 
A. DPD 

DPD was introduced in 1992 by Hoogerbrugge and Koel- 
man lfl5l [l6ll as a particle based simulation technique for 
complex fluids. It was originally considered as a blend be- 
tween the Molecular Dynamics and Lattice Gas Automata 
techniques, the first because of the similar integration scheme 
and representation, and the latter because it is focused mainly 
on respecting the mechanical conservation laws. 

The DPD equations of motion, with particle positions r,, 
velocities v, and momenta p,, can be written as a system of 
Langevin equations 
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(1) 
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where F?, F? and F? are conservative, dissipative, and 
stochastic forces between particles i and j. A particular prop- 
erty of DPD is that all interactions are modeled by pairwise 
central forces. This ensures the conservation of linear and 
angular momentum which in turn implies the correct Navier- 
Stokes equations in the continuum limit lfl7l[l8ll . 

In the original formulation of DPD, the dissipative and 
stochastic forces were chosen independently. In what later 
has become the standard DPD model, formulated by Espanol 
and Warren lfl9ll . the dissipative and stochastic forces are not 
independent; rather, they depend on each other through a fluc- 
tuation dissipation theorem and together function as a local 



pairwise thermostat. The equilibrium properties of a DPD sys- 
tem are determined by the conservative force alone, while the 
thermostat defines how the system reaches equilibrium. The 
common way to express the dissipative and stochastic forces 
is 
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where ry is the distance between particles ; and j, ey is the 
unit vector pointing from j to i, and Vy is the velocity dif- 
ference Vy = v j — Vj. £b is Boltzmann's constant, and T is 
the temperature. The scalar function ffl(ry) describes how 
the stochastic and dissipative forces depend on the distance 
between the particles, and Qj is interpreted as a symmetric 
Gaussian white noise term with mean zero and covariance 



(5) 



where 5y and 5(f) are the Kronecker and Dirac delta func- 
tions, respectively. 

It has become standard practice to choose the conservative 
and stochastic interactions in DPD to be linearly decreasing 
functions with a finite cutoff radius. To give the simulated 
system a physical interpretation the magnitude of the conser- 
vative force is generally tuned such that some desired equi- 
librium property of the system, e.g. the compressibility 12011 . 
is matched. The magnitude of the stochastic force is then set 
such that the temperature equilibration is reasonably fast, but 
still slow enough to allow a stable simulation. 

We and others have previously argued that the DPD ansatz 
is to be seen as the result of a systematic coarse-graining from 
a given microscopic system and that neither the conservative 
nor the dissipative/stochastic interactions can be chosen arbi- 
trarily JH, IIH H2]. Rather, the DPD interactions should re- 
sult from projecting the microscopic dynamics onto a coarse- 
grained level. This argument is valid regardless of the pro- 
jection, as long as the projected dynamics truly follows the 
DPD ansatz. Types of projections that can be considered are 
for example clustering collections of freely moving particles 
into sin gle DPD beads, which is the way DPD is commonly 
viewed |H [M IH HI HI |H, or DPD used as a thermostat 
for united atoms coarse-graining JH H^, l27ll . In the former 
case, the stochastic and dissipative interactions result from the 
internal motion of the microscopic particles inside the DPD 
beads, whereas for the latter case, these interactions result 
from the internal atomistic motion not explicitly simulated in 
the coarse-grained particles. 



B. Deriving the DPD forces 

The effective interactions in the DPD model is derived us- 
ing the trajectories from a detailed molecular dynamics sim- 
ulation. The trajectories of the particles in the united atoms 
representation are calculated as a center of mass projection 
from the atomistic trajectories. The goal of our method is to 
find forces acting on the center of mass point particles that 
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give equilibrium distribution and transport properties that are 
as close as possible to those caused by the true force field. 

We use a tiered strategy, where we first find the radial distri- 
bution function (RDF) of the united atoms system, calculated 
from the simulations of the underlying SPC system. Given 
this function, we use the Inverse Monte Carlo (IMC) method 
JH 0] to find the potential function that gives rise to the same 
RDF. This and similar methods rest on the result by Hender- 
son 0], who showed the uniqueness of such a potential. The 
IMC method is an iterative method which may be viewed as 
a multi-dimensional Newton-Raphson technique. In each it- 
eration the potential is updated using the gradient of the RDF 
with respect to the potential, which is estimated from standard 
forward Monte Carlo simulations using the present potential. 

Once the IMC method has converged to a solution, the con- 
servative force is fixed. In order to complete the DPD descrip- 
tion of the united atoms system, it remains to determine the 
dissipative and stochastic forces. Given that these two forces 
obey a dissipation-fluctuation relation, the choice of stochastic 
force does not influence the equilibrium distribution, only the 
approach to equilibrium and hence, in general, the transport 
properties of the system. Unfortunately, this also means that 
the RDF gives us no information about the relevant stochas- 
tic and dissipative forces, and we must turn to observables 
that more directly reflects how the shape and magnitude of 
the stochastic force determines the transport properties. 

In lfl3ll we have suggested a practical scheme based on the 
auto-correlation of the force fluctuations, K F (r, t) , defined as 



K F (r,T) = (8,(S F pi)-8,(8 F pj)) r , 



(6) 



where (} r denotes ensemble average conditional on the dis- 
tance r between particles ; and j at time zero, and where 

8 t (8 F p,) = p,(f) -p,(0) - fdt [if (»)+!?(»)] (7) 



is the contribution to particle i s change in momentum from 
the stochastic forces during the time interval [0,t]. [In our 
notation 8 t denotes a difference in time, and 8 F the difference 
between total force and deterministic (conservative and dissi- 
pative) forces.] 

Consider now the thought experiment of adding a force to 
the SPC simulation to freeze the center of mass of each wa- 
ter molecule in place at time zero, but allowing the internal 
degrees of freedom of each molecule to evolve as usual (i.e. 
the orientation and internal vibration of the molecules). Us- 
ing Mori -Zwanzig theory one can then show the exact relation 



(0 2 (r) = lim - 
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2k B T dx 



K F (r,r) 
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where ffl(r) is the function determining both the dissipative 
and the stochastic DPD interactions. In the SPC simulations, 
no such force is added, but Eq. (O may still hold approx- 
imately in a region of T if the internal degrees of freedom 
change on a faster time scale than the center of mass positions 
of the molecules, i.e., if there is a region where K F (r, t) is an 
approximately linear function of time. In lfl3ll we showed that 



this is indeed the case if the DPD ansatz holds. In the underly- 
ing SPC system, the DPD assumption that the time autocorre- 
lation of the center of mass forces is negligible does not hold 
for short time scales. In this case, the time region must be 
chosen sufficiently long that the force has time to decorrelate, 
but still short enough that the particles do not move too far 
from their positions at the beginning of the time interval. 

Since the right-hand side of Eq. ([8]) depends on the func- 
tion 0)(r), our method is iterative. As a starting point, we 
take co(r) = 0, corresponding to a strictly deterministic time 
evolution of the center of masses. The method is then as fol- 
lows: First, for each value of r evaluate K F (r, t) as a function 
of time. Second, identify the time region where K F (r, t) is 
approximately linear (this region may be different for differ- 
ent values of r), and estimate co(r) as the right-hand side of 
Eq. ([8]l. If the new and the old CO(r) are close for all values of 
r, Eq. ([8]) holds approximately and we may stop. Otherwise, 
we repeat the procedure from step one. This concludes our 
background section. Before we present our results, in the next 
section we give the details of the molecular dynamics simula- 
tion of the SPC model. 



III. SIMULATION SETUP 

We use Simple Point Charge (SPC) water as the model sys- 
tem. The SPC model is not the most accurate water model 
available when it comes to representing transport proper- 
ties JUt], but is widely used and serves as a good test case for 
our method. The simulations were performed with a periodic 
cubic box of size 3.0 nm, with a total of 895 SPC molecules, 
using the molecular dynamics software GROMACS. The inte- 
gration time step was set to 2.0 fs, and temperature and pres- 
sure equilibration was performed using the Berendsen ther- 
mostat, with target temperature and pressure 298 K and 1 
bar respectively. The compressibility parameter was set to 
4.6 x 10~ 5 (1/bar). Long range electrostatic forces was han- 
dled using the reaction-field method, with a cut-off of 1.4 
nm. From an initial simulation, the radial distribution func- 
tion (RDF) of the center of mass of the SPC molecules was 
measured. This data was then used as input to the inverse 
Monte Carlo method, which produces an effective potential 
for the center of mass motion of the coarse-grained represen- 
tation of the SPC system. Given this potential the GROMACS 
SPC simulation was reran, now measuring k f [cf. Eq. (0], 
and gradually changing co(r) until convergence. 



IV. RESULTS 
A. Conservative potential 

Following the procedure in we took the potential of 
mean force, 4>mf = — ksTlog[g(r)], as the initial guess for the 
effective potential, shown as the dashed graph in Fig.Q] The 
final effective potential obtained using IMC, is shown as the 
solid graph in the same figure. This is the result after 10 itera- 
tions, each using 5 million Monte Carlo updates. The resulting 
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FIG. 1 : Pairwise potential derived from the radial distribution func- 
tion for the center of mass of the SPC molecules by using the IMC 
method (solid). The potential converges to zero at about r = 1.0 nm. 
The potential of mean force is plotted as a reference (dashed line), 
deviating substantially from the estimated potential. 

potential is clearly different from the potential of mean force, 
and has the important property that it faithfully represents the 
equilibrium properties of the system according to the theorem 
of Henderson |4[]. This is true in any simulation, independent 
of the thermostat, as long as the thermostat is constructed ac- 
cording to a fluctuation-dissipation theorem Il3ll . 

B. Dissipative/Stochastic forces 

The potential derived in the previous section determines the 
functional form of the conservative forces, Ffj, between the 
particles. In a simulation, the effective conservative force on 
a coarse-grained particle can then be calculated from Ff = 
Hj^iFi]( r ij)- From Eq. ^ we get an estimate of $(5f-p,-) 
under the initial assumption that the dissipative forces are all 
equal to zero. By using Eqs. (H)) and ([6), a new estimate of 
Co(r) can be found. Inserting this new approximation in the 
dissipative force in Eq. @, and repeating the procedure until 
convergence produces the co(r) function that should be used 
together with the conservative potential in a DPD representa- 
tion of the original SPC system. 

Fig. (O shows the measured values of Kp for the first itera- 
tion step. For small times, the motion of the SPC molecules' 
mass centers are smooth and deterministic as determined by 
the forces in the detailed SPC model. On this time scale, the 
white noise approximation of the DPD method is ill-defined. 
It is only on longer time scales that the degrees of freedom 
taken away by the projection can be assumed to affect the 
slow (center of mass) degrees of freedom in the form of noise 
and dissipation. At these time scales, we expect from the- 
ory to find a linear region in Kp lfl3ll . and estimate (o(r) from 
the slope [cf. Eq. ((§)]. We identify this as the region above 
T = 0.3 ps. For times significantly larger than T = 0.6 ps, the 
center of mass of the SPC molecules will have had time to 
move a non-negligible distance from their original positions, 
rendering the conditioning on r in Eq. © ill-defined. The re- 
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FIG. 2: The asymptotic slope of Kf defines, through Eq. (8), the 
functional form of £i)(r). In this figure, the values of Kp are shown 
as a function of the time difference, T [cf. Eq. 0], for four different 
distance values: r = 0.27 nm (x), r = 0.33 nm (□), r = 0.39 nm (o), 
and r = 0.45 nm (o). The solid lines are linear fits of the Kp function 
in the range T £ [0.3, 0.6] ps. 



gion from which the slope of Kp can be measured must hence 
be chosen smaller than this time. In our simulations the SPC 
molecules move on average 0.14 nm during 0.6 ps, which is 
close to half the distance to the first peak in the RDF. Hence, 
we conclude that the particles stay close to their initial posi- 
tion during the time interval used to determine ffl(r). 

The resulting form of <a(r) is shown in Fig. [3] The curves 
represent the first and last steps in the iterative procedure. Af- 
ter eight updates of ffl(r), the procedure has converged to the 
function shown as the solid line. The inset shows the norm of 
the distance between co(r) of consecutive iterations, defined 
as 

||Ao%|| = f (flfc(0-^-i(0)Vj£^(OL<i(0> (9) 

i=0 V i=0 i=0 

where k is the iteration number, and ; goes over all points on 
the co-curves. 

Using the derived interactions, i.e., the conservative poten- 
tial and CO(r), in a DPD model, we can measure the transport 
properties of the coarse-grained model. Table U presents the 
resulting diffusion and viscosity for the DPD model, com- 
pared to the original SPC simulation and also to the simple 
united atoms (UA) model with no dissipative and stochastic 
forces. We note that both the diffusion and viscosity are off 
by roughly a factor of four in the united atoms simulation, 
whereas the addition of the newly derived DPD thermostat re- 
duces that number to 60 percent for the diffusion and only 20 
percent for the viscosity. 

In Fig.@]we compare the newly derived co(r) function from 
the bottom-up method presented in this paper with the top- 
down approach in reference ifTill . The dashed curves repre- 
sent hand tuned (o(r) functions that all give approximately 
correct diffusion and viscosity, compared to that of the de- 
tailed SPC simulation. The solid line is the co(r) function 
derived without any manual tuning in this paper. The rather 
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FIG. 3: The figure shows ffl(r) as a function of the inter-particle dis- 
tance r, after the first (dashed line) and last (solid line) steps of the 
iterative procedure. Starting with ffl(r) = for all r, the slope of Kf 
is measured (cf. Fig.O, and a new estimate of co(r) is obtained by 
using Eq. l[8j. The convergence towards a self consistent solution is 
shown as inset, with Aa)^ defined in Eq. ((9). After eight steps, a self 
consistent solution is reached, resulting in a Aco-value of approxi- 
mately 0. 

Method D [10~ 9 m 2 s' 1 ] 77 [cP] 
SPC 3.87 ± 0.04 0.44 ±0.05 
UA 15.20±0.14 0.120±0.005 

DPD 6.13 ±0.07 0.36 ±0.05 

TABLE I: This table presents a comparison between the diffusion 
coefficient and the viscosity for the microscopic system (SPC) and 
for two coarse-grained representations of the same system (UA and 
DPD) described in the main text. A clear improvement of the trans- 
port properties of the coarse-grained system can be seen when adding 
the DPD thermostat. 



close agreement between transport properties of the coarse- 
grained (DPD) model and the SPC model presented in table U 
is further strengthened by this figure, as we note that the dif- 
ference between the derived Co(r) function and the hand tuned 
curves is relatively small, though it is worth noting that the 
derived DPD interactions is consistently weaker than the hand 
tuned thermostat. 



V. DISCUSSION 

We conclude with a brief summary and discussion of our 
results. First, we have calculated the effective conservative 
force for a united atoms representation of water from the ra- 
dial distribution function of the center of mass. Simulations 
using the resulting forces show that the transport properties of 
the united atoms model differ significantly from those of SPC 
water: the diffusion rate is about four times too high, and the 
viscosity is approximately four times too low. This discrep- 
ancy is not too surprising considering that we have replaced 
the strongly polar SPC water molecule with a simple point 
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FIG. 4: The figure shows a comparison between the co(r) function 
derived in this paper and some examples of the functional forms of 
(0(f) that were manually tuned to give approximately correct diffu- 
sion and viscosity in reference fl4Tl . 



particle with a radially symmetric interaction potential. 

Second, our iterative method for estimating the stochastic 
force in the DPD extension of the united atoms model finds 
the force that best fulfills the self-consistency relation between 
Kf and co. Including the resulting dissipative and stochastic 
forces gives transport properties that are significantly closer 
to the SPC simulations. For instance, the DPD diffusion rate 
is approximately 60% higher than that of the SPC simulations, 
but excluding the dissipative and stochastic forces yields a dif- 
fusion rate that is 290% higher (cf. Table H). To give some 
perspective to these figures, consider that the difference be- 
tween the experimental value of the diffusion rate of water 
(2.27 x lO- 9 m 2 /s JH) and the diffusion rate of SPC is of 
the same order as the difference between the diffusion rates of 
SPC and DPD. 

Using a linear ansatz for 0)(r), we showed in a previous 
report 1 14-] that it is possible to obtain a simultaneous agree- 
ment of the diffusion rate and viscosity by carefully choosing 
the support and magnitude of £»(r). Since these are defined 
on macroscopic time and space scales, it follows that match- 
ing these properties does not necessarily preserve the dynam- 
ics at the mesoscale. It is symptomatic that several different 
parameter settings was found that gave similar transport prop- 
erties 11411 , The transport properties using the a>(r) found in 
this article are not as accurate, but on the other hand no tuning 
was made of 0)(r) to improve the transport properties. The 
method presented in this paper is entirely a "bottom-up" ap- 
proach, as opposed to the "top-down" method used in stag- 
nated DPD and in fl4Tl . Given the constrained nature of the 
center of mass representation of water, it is not unexpected 
that we may have to choose which transport property to best 
represent. The advantage of the method in this paper is that it 
does not make any assumptions about the shape of (0{r). If the 
DPD ansatz is approximately valid, the method automatically 
identifies the shape, magnitude and support for the stochastic 
force, without the need to tune any parameters. 
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The water molecule presents a difficult case for united 
atoms modeling. Because the dipole-dipole interactions are 
strong and fluctuate on relatively long time scales (the dipole- 
dipole decorrelation time is approximately 4.7 ps ll30ll . much 
longer than the time over which we estimate the stochastic 
forces), the center of mass motion of SPC water is neither av- 
eraging nor a Langevin process. This suggest that we have 
underestimated the role of the fluctuations. However, the es- 
timated conservative interaction includes information of typ- 
ical local configurations, and the stochastic interactions cap- 
ture fluctuations around the effective mean forces. The fluc- 
tuations may therefore very well equilibrate around the local 
mean configurations much faster than the total decorrelation 
time between the local orientation of the molecules. We think 
that this is the explanation for the linear region of Kf to appear 



on relatively short time scales. 

Furthermore, we have seen that including the stochastic 
forces gives a significant improvement over only the conser- 
vative forces. This shows that it is possible to find DPD forces 
that capture not only the equilibrium distribution but also the 
transport properties of fluids with a rich spatial structure and 
time evolution, without explicitly representing the polar dy- 
namics. As discussed in the introduction, this can be useful 
for simulating water in the bulk, where the main concern is 
the water-water interactions. 
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